# working directory ---------------------------
setwd("~/Library/CloudStorage/OneDrive-UniversitédeGenève/oecd-dac-project-unige/paper-gender-targeting/2025-ddfi-isq-replication")

# packages 
library(tidyverse)
library(ggeffects)
library(reshape2) # for melt() and dcast()
library(gridExtra)

# figure 1 ---------------------------
f1_data <- read_csv("2025-ddfi-isq-data-f1.csv")

ggplot(f1_data, aes(x = y, y = gepm_pct, fill = gepm_cat)) + 
  geom_col() + 
  geom_point(aes(y = commit_purpose_both_pct)) + 
  theme_minimal() + 
  theme(aspect.ratio = 5/16, 
        text = element_text(size = 18), 
        axis.title.x = element_text(vjust = -1),
        axis.title.y = element_text(size = 15),
        legend.position = "top", 
        legend.direction = "vertical") +
  scale_x_continuous(n.breaks = 15) +
  scale_y_continuous(labels = scales::percent) +
  scale_fill_manual(values = c("lightgray", "darkgray"), labels = c("Gender Mainstreaming", "Gender Equality as Primary Target"), name = "New Legend Title") +
  guides(fill = guide_legend("Black dots indicate % of ODA commitments to gender-related sectors")) + 
  xlab("") +
  ylab("% ODA Commitments From All DAC Donors") -> f1_a

f1_data %>% 
  filter(gepm_cat == "commit_gender_2_pct") %>% 
  ggplot(aes(x = y, y = gepm_pct, fill = gepm_cat)) + 
  geom_col() + 
  geom_point(aes(y = commit_purpose_both_pct)) + 
  theme_minimal() + 
  theme(aspect.ratio = 5/16, 
        text = element_text(size = 18), 
        axis.title.x = element_text(vjust = -1),
        axis.title.y = element_text(size = 15),
        legend.position = "none", 
        legend.title = element_blank()) +  
  scale_x_continuous(n.breaks = 15) +
  scale_y_continuous(labels = scales::percent) +
  scale_fill_manual(values = "darkgray", name = "New Legend Title") +
  xlab("") +
  ylab("% ODA Commitments From All DAC Donors") -> f1_b

grid.arrange(f1_a, f1_b, ncol = 1, heights = c(1.24, 1)) -> f1
ggsave("2025-ddfi-isq-f1.pdf", plot = f1, height = 10, width = 16, units = "in")


# figure 2 ---------------------------
f2_data <- read_csv("2025-ddfi-isq-data-f2.csv")

f2_data %>% 
  rename(Bypass = commit_bypass_pct, Public = commit_public_pct, "NA" = commit_na_pct) %>% 
  melt() %>% 
  ggplot(aes(x = type, y = value, fill = variable, label = scales::percent(value, accuracy = 1))) + 
  geom_col(width = .75) + 
  geom_text(position = position_stack(), vjust = -.35, size = 5, color = "#3b3b3b") + 
  theme_minimal() + 
  theme(aspect.ratio = 10/16, 
        text = element_text(size = 18), 
        legend.position = "top", 
        legend.title = element_blank(), 
        legend.direction = "horizontal") +
  scale_y_continuous(labels = scales::percent) +
  scale_fill_manual(values = c("lightgray", "darkgray", "#3b3b3b")) +
  xlab("") +
  ylab("% ODA Commitments From All DAC Donors") -> f2
ggsave("2025-ddfi-isq-f2.pdf", plot = f2, height = 10, width = 16, units = "in")


# figure 3 ---------------------------
ry <- read_csv("2025-ddfi-isq-data-ry.csv") # recipient-year level, built including eu lines 

ry %>% 
  mutate(gdppc_log_ry = log(gdppc_ry), 
         quota_adpt_ry = if_else(quota_adpt_ry == 0, "No", "Yes"),
         vdem_row_id_ry = if_else(vdem_row_id_ry == "democracy", "Democracy", "Autocracy")) -> ry

m3 <- lm(commit_gender_1_pct ~ vdem_row_id_ry * quota_adpt_ry
           + gdppc_log_ry + dependence_ry + y,
           data = ry)

ggeffect(m3, terms = c("quota_adpt_ry", "vdem_row_id_ry"), 
         vcov_fun = "vcovCL", vcov_type = "HC1",
         vcov_args = list(cluster = ry$r),
         robust = TRUE) -> data_f3_m3

ggplot(data_f3_m3, aes(x, predicted)) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high, colour = group), position = position_dodge(width = 0.25)) +
  scale_color_manual(values = c("black", "darkgray")) +
  theme_minimal() +
  theme(aspect.ratio = 16/16, 
        text = element_text(size = 20), 
        legend.position = "top", 
        legend.title = element_blank(), 
        legend.text = element_text(size = 20),
        legend.direction = "horizontal") + 
  xlab("Quota") + 
  ylab("% Mainstreamed Commitments") -> f3a
  
m4 <- lm(commit_gender_1_pct ~ vdem_row_id_ry * wbl_ry 
         + gdppc_log_ry + dependence_ry + y, 
         data = ry)

ggeffect(m4, terms = c("wbl_ry", "vdem_row_id_ry"), 
         vcov_fun = "vcovCL", vcov_type = "HC1",
         vcov_args = list(cluster = ry$r),
         robust = TRUE) -> data_f3_m4

ggplot(data_f3_m4, aes(x, predicted)) +
  geom_line(aes(linetype = group)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.15) +
  scale_linetype_manual(values = c("solid", "dashed", "dotted")) +
  scale_fill_manual(values = c("black", "darkgray")) +
  theme_minimal() +
  theme(aspect.ratio = 16/16, 
        text = element_text(size = 20), 
        legend.position = "top", 
        legend.title = element_blank(), 
        legend.text = element_text(size = 20),
        legend.direction = "horizontal") + 
  xlab("WBL Index") + 
  ylab("% Mainstreamed Commitments") -> f3b

grid.arrange(f3a, f3b, ncol = 2) -> f3
ggsave("2025-ddfi-isq-f3.pdf", plot = f3, height = 10, width = 16, units = "in")


# figure a1 ---------------------------
ry %>% 
  filter(y %in% c(2018:2022)) %>% 
  mutate(r = if_else(r == "China (People's Republic of)", "China", 
                     if_else(r == "Syrian Arab Republic", "Syria",
                             if_else(r == "Viet Nam", "Vietnam", r)))) -> ry_5yrs
recipients <- sort(unique(ry_5yrs$r))

# effective sample:
# what recipients stay in the data for models 3-4? 
ry_5yrs %>% 
  filter(!is.na(commit_gender_1_pct), 
         !is.na(vdem_row_id_ry),
         !is.na(quota_adpt_ry),
         !is.na(wbl_ry),
         !is.na(gdppc_ry),
         !is.na(dependence_ry)) -> ry_5yrs_effective 
recipients_effective <- sort(unique(ry_5yrs_effective$r))

setdiff(recipients, recipients_effective)
paste(setdiff(recipients, recipients_effective), collapse = ', ') 

# on average, who got the most and least GEPM aid over the past five years? 
# only recipients in effective sample 
ry_5yrs %>% 
  mutate(in_sample = if_else(r %in% recipients_effective, 1, 0)) %>% 
  filter(in_sample == 1) %>% 
  group_by(r) %>% 
  summarize(commit_gender_1_pct_avg = mean(commit_gender_1_pct),
            commit_gender_2_pct_avg = mean(commit_gender_2_pct), 
            commit_gender_both_pct_avg = mean(commit_gender_both_pct), 
            commit_purpose_both_pct_avg = mean(commit_purpose_both_pct)) %>% 
  ungroup() %>% 
  mutate(commit_gender_1_pct_avg_rank = rank(commit_gender_1_pct_avg),
         commit_gender_2_pct_avg_rank = rank(commit_gender_2_pct_avg),
         commit_gender_both_pct_avg_rank = rank(commit_gender_both_pct_avg),
         commit_purpose_both_pct_avg_rank = rank(commit_purpose_both_pct_avg)) -> ry_5yrs_avg 

# all recipients, ranked by GEPM 1 + 2 
top_bottom_15 <- c(1:15, 102:116)
ry_5yrs_avg %>% 
  filter(commit_gender_both_pct_avg_rank %in% top_bottom_15) %>% 
  select(r, commit_gender_both_pct_avg_rank, 
         commit_purpose_both_pct_avg, commit_gender_1_pct_avg, commit_gender_2_pct_avg) %>% 
  melt(id.vars = c("r", "commit_gender_both_pct_avg_rank","commit_purpose_both_pct_avg")) %>% 
  ggplot(aes(x = value, y = reorder(r, commit_gender_both_pct_avg_rank), fill = variable)) + 
  geom_col() + 
  geom_point(aes(x = commit_purpose_both_pct_avg)) +
  theme_minimal() + 
  theme(aspect.ratio = 10/10, 
        text = element_text(size = 18), 
        legend.position = "top", 
        legend.direction = "vertical") +
  scale_x_continuous(labels = function(value) paste0(value, "%")) +
  scale_fill_manual(values = c("lightgray", "darkgray"), labels = c("Gender Mainstreaming", "Gender Equality as Primary Target")) +
  guides(fill = guide_legend("Black dots indicate % of ODA commitments to gender-related sectors"), size = 10) + 
  ylab("") +
  xlab("% ODA Commitments From All DAC Donors") -> fa1
ggsave("2025-ddfi-isq-fa1.pdf", plot = fa1, height = 10, width = 10, units = "in")